Spatial decay of the single-particle density matrix in insulators: 
analytic results in two and three dimensions 
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Analytic results for the asymptotic decay of the electron density matrix in insulators have been 
obtained in all three dimensions {D = 1 — 3) for a tight-binding model defined on a simple cubic 
lattice. The anisotropic decay length is shown to be dependent on the energy parameters of the 
model. The existence of the power-law prefactor, oc r"^''^, is demonstrated. 

PACS numbers: 71.15.Ap, 71.20.-b, 71.15.-m 



In a phenomenological classical approach to atomic dy- 
namics, a very local picture of interatomic interactions is 
often highly appropriate. The simplest example is the 
use of Keating "springs" to describe small atomic oscil- 
lations in solids [Q. In quantum mechanics (QM), the 
situation is superficially quite different. If are single- 
particle wave functions (for example, the Kohn-Sham or- 
bitals of density-functional theory), then the electronic 
contributions to the total energy, forces or the dynamical 
matrix can all be obtained from the single-particle den- 
sity matrix (DM): p(r,r') = Eocc V'l* (r)V'i(r') §. For 
a condensed-matter system, virtually all of the Hamilto- 
nian eigenstates ipi are extended and oscillatory through- 
out the volume of the system, except some of the states 
near the band edges in the case of disordered systems. 
To the extent that p(r,r') is built from objects that are 
extended in real space, it is a nontrivial feature of the 
quantum mechanics of extended systems that in fact the 
DM can be localized, even exponentially so in insulators 

ii- 

The possibility of a local formulation of QM goes back 
at least to Wannier (5 , whose " Wannier functions" decay 
in real space. Kohn [ 3| extended this work, and empha- 
sized the principle of " nearsightedness" , that is the de- 
pendence of local properties, such as forces, on the local 
environment. The possibility of real-space local formu- 
lations of QM has inspired vigorous activity on efficient 
"order-N" methods for electronic structure 

In early pioneering work, Kohn |Q showed that the DM 
and Wannier functions for a ID centrosymmctric two- 
band model decay exponentially in insulators: p(r, 0) cx 
exp(— A|r|), where A cx E^ip 0, although Ismail-Beigi 
and Arias 1^ find instead that A cx Eg^p as Eg^p 0. Dcs 
Cloizeaux proved that the DM decays exponentially quite 
generally in insulators [^. Recently, He and Vanderbilt 
pof demonstrated that in ID there is in fact a power- 
law prefactor cx r~^/^, and that the characteristic inverse 
decay length. A, of both the DM and Wannier functions 
is universal, although the prefactor is different in each 
case. 

In this Letter, we report analytic results for the spa- 



tial decay of the DM for a two-band tight-binding (TB) 
Hamiltonian for all three dimensions, D — 1 — 3. The 
Hamiltonian we use is similar in spirit to that of Harri- 
son's bond orbital model [0 and is the minimal model 
which contains the basic features of an insulator. Both 
the power-law prefactor and the exponential spatial de- 
cay of the DM are found in the asymptotic regime. Thus, 
we reproduce the results of Refs. [0, |o| , and extend them 
to 2D and 3D. We derive the dependence of the inverse 
decay length. A, in terms of the fundamental energy pa- 
rameters of the model. We further show how, in certain 
cases, the exponential localization arises from the trun- 
cation of the sum defining p. Precise numerical results 
support the analytic asymptotic behaviour in all three 
dimensions. 

Our model is of centrosymmctric nearest-neighbor TB 
form on a simple cubic (s.c.) lattice described by the 
Hamiltonian: 



(1) 



Here, each site i is characterized by two bare orthogonal 
electronic states, (/i = 1,2) with bare energies 
(Ae = El — £2 > for definiteness) , which together form 
the orthonormal bare basis set. These may be interpreted 
as bonding and antibonding states centered at each lat- 
tice site. The hopping integrals, t^^, between similar 
states (the same n) on different nearest-neighbour sites 
j (to i) transform the on-site energy levels into bands, 
while the hopping integrals, t^^i , between different states 
(/i 7^ ^') on different sites are responsible for inter-band 
hybridization. For this model, the eigenstates |k,7) (k 
is the wavevector and 7 = 1,2 indexes the bands) and 
dispersion are analytically available and the spectrum 
exhibits the main features of a semiconductor /insulator 
(i.e. two bands separated by a gap for a certain range of 
parameters). 

We have written the density operator for the Hamil- 
tonian given by Eq. (|l|), first, in the eigenvector basis, 
|k, 7), then expanded the eigenstates in the site basis 
\ifi) and evaluated the density matrix elements, Pifi.j^' 
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FIG. 1: The dependence on distance of the logarithm of the 
DM llnp^^l (a,b) and of the effective inverse decay length, 
A — — ln|p^^|/cr (c,d) in the 2D case (the prefactor is in- 
cluded in the definition of A for easier comparison of analyti- 
cal results with numerical ones). In (a) and (c), the DM has 
been calculated on the lattice sites situated around the main 
diagonal [11] and characterized by polar angles 4> € [lO^jSO"], 
while in (b) and (d), Vy — and i-'x varies (in the [10] direc- 
tion). The circles represent exact numerical results obtained 
using Eq. ^ or Eq. (^), while the crosses are the results of 
the analytical expressions (§)-(^ (they are practically indis- 
tinguishable in (a) and (c)). In (c), the different analytical 
contributions to the effective inverse decay length are shown 
by the solid line (the ln(2yl)-term in Eq. (^)), squares (the 
second term in Eq. (^), triangles (the sum of the last two 
terms in Eq. (^) and diamonds (the prefactor in Eq. (^). In 
(d), the solid line is the ln(2AV-contribution to the effective 
inverse decay length (see Eq. (|6|)), the squares represent the 
contribution of (— In | J^^ (z/x/LA)|/j^x) to A and the diamonds 
signify the prefactor in Eq. 



in the site basis. For fixed site indices i and j, the DM 
contains four elements, which, in the case of a semicon- 
ductor at zero temperature, obey the following relations: 
Pit^,m = ^ii s-'^d piij2 = Pi2ji- Below we investigate 
the decay properties only for the off-diagonal elements, 
Pii,j2 (the analysis and results for the diagonal elements 
Pifj.,jfi are very similar to those for Piij2 presented below). 
The expression for pa,j2 = p{^ij) can be conveniently 
represented via an integral in reciprocal space over the 
first Brillouin zone: 



'5k 



(2) 



where {ai^a;, avy^ avz} — — is the connection vec- 
tor (with the unit-cell parameter a taken hereafter to be 
unity, and (e.g. a = x,y, z m 3D) being an inte- 
ger), S'k = (1/2) ^^-^-^ exp (ik • Tij) is the structure fac- 
tor and = (Ae -f 2(tii — t22)Sk)/Ati2 is the only en- 
ergy parameter of the model which, for the symmetric 
case (^11 = <22 = t), analyzed below for simplicity, is 



k- independent i.e. Ak = A = Ae/4ti2. In the model 
without interband hybridization (ti2 = or A — > cxd), 
the DM given by Eq. (||) is obviously pij = 0. In order to 
study the spatial decay of the density matrix, we should 
take the integral in Eq. ^ and investigate its depen- 
dence on the lattice indices i^a ■ This can be done analyt- 
ically in the asymptotic regime for all or some i/q, 1, if 
the energy parameter of the model obeys the inequality, 
A > D, which is the case for weak interband hybridiza- 
tion in semiconductors. 

Eq. (|^) has been derived without making use of the 
particular symmetry type of the underlying lattice, and 
is valid for any primitive lattice (one atom per unit cell) . 
The structure factor, S'k, reflects the differences between 
various lattices. The DM given by Eq. (^) can be eas- 
ily calculated numerically for any reasonable choice of 
parameters, and the results of such calculations support 
the main conclusions derived below analytically for the 
s.c. lattice, subject to some additional restrictions (sym- 
metric case, til = ^22, and weak hybridization, A > D). 

To proceed with the evaluation of the DM given by 
Eq. (^, first we expand in a series the denominator of 
the integrand (S^/A < 1, if A > D), separate the vari- 
ables and then evaluate all D integrals making use of the 
orthogonality of the cos(fcct;/Q.)-functions involved. This 
results in the following general exact form for the DM for 
all three dimensions: 



P^c 



(-ir 



(4/1)2^- 



-1 /l^ 

fe=0 



(2fc')' 



{2k' + i)j:d 



h^y.{k - m + Vy)!] be- 
but that for 31?, 
to)! (to -|- fx)!(m + 



(4:A)''{k' 

(3) 

if <T = h'a is odd, and zero otherwise. Here k' ~ 
k + V = k + {a ~ l)/2. El = l/[fc!(fc + v^y.] and 
^2,3 are hypergeometric finite sums, the one for 21?, 
^2'= ELol/H(fc-m)!(m 
ing expressible in closed form |lj 
^3^j:t=oi'^rn + u^ + i.yy./[mlik 
t'y)!(TO + I'-x + ^'y)!(^ — m + J^z)!], cannot so be written, 
in principle |l^ . The orthogonality of the cos(fcct^'a)- 
functions, being a consequence of the Fermi level lying 
in the gap (only in this case are the upper limits of the 
integrals in Eq. (H) equal to tt), is an important feature 
in deriving Eq. (0). The original series is truncated for 
small indices which results in the appearance of the ex- 
ponential factor, (4A)~2''-i^ before the sum in Eq. (|^). 
This factor then appreciably contributes to the energy 
(A)-dcpendent part of the decay length of the DM (al- 
though see Fig. |^(e) where this is not the case). 

Further analysis is based on the use of Stirling's ap- 
proximation for Va ^ 1 which allows us to evaluate 
asymptotically both the internal finite sum S3 in the 3-D 
case and the external infinite sums over k in all dimen- 
sions. For example, if all ^ 1, then keeping the lead- 
ing term in the asymptotic expansions, we can write the 
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FIG. 2: The dependence on distance of | In p^^ \ (a-c) and 
of A = — ln|p^^|/(T (d-f) in the 3-D case. In (a) and (d), the 
DM has been calculated on the lattice sites situated around 
the main diagonal [111] for solid angles defined by the polar 
angle (j) G [20^ 70°] and azimuthal angle 6 £ [20°, 70% in (b) 
and (e), ly^ — and and i^y vary around the plane diag- 
onal and are characterized by the polar angle (f) G [20°, 70°] 
(around the [110] direction); and in (c) and (f), = fy = 
and varies (in the [100] direction). All the symbols have 
the same meaning as in Fig. |l| except in (e), where the solid 
circles represent the contribution to A due to the exponential 
term in Eq. (^, the squares and triangles are due to J^^ and 
"^fx+fyi respectively, while the diamonds signify the prefactor 
in Eq. d). 



following expression for the DM in all three dimensions: 

p.„^(-l)''(27r)-^/Mn^.J exp[-aA] , (4) 

where the inverse decay length, X{A, Ua), is given by 

A(A^^a) = ln(2A) + ^ + ^-lna , (5) 



with cr_i = Y.'^a^ and CTin = Va\nva. As fol- 
lows from Eqs. (0)-(||) and Figs. 0(a,c) and §(a,d), the 
DM in this model of an insulator: (i) exponentially de- 
cays with the distance-related parameter cr; (ii) has a 
power-law prefactor which is proportional to v~^l'^\ (iii) 
is anisotropic; and (iv) the inverse decay length of the 
DM depends on the energy parameters of the model. The 
behaviour with effective distance of In \pu^ \ and the effec- 
tive inverse decay length, A = — In l/ci, are shown in 
Figs. |l|(a,c) for the 2D case and in Figs. ||(a,d) for the 
31? case, respectively. The different points correspond to 
all the lattice nodes within a relatively wide solid (po- 
lar in 2D) angle around the main diagonal. Examina- 
tion of the results shows that the slowest decay of the 
DM occurs along the main diagonal, i.e. the [111] direc- 
tion in the 313 case. The analytical expression ^ allows 
the different contributions to the effective decay length 




FIG. 3: The dependence of the inverse decay length A (a) 
and the effective inverse decay length A (b) on the energy- 
dependent parameter A of the model in 3-D. The data in (a) 
given by the thick solid and long dashed lines were obtained by 
linear regression from the slope of In |p| evaluated numerically 
using Eqs. (^ and (^ for the [111] direction, respectively. The 
thin solid line in (a) represents Eq. (^) (the difference between 
the thin and thick solid lines is due to the prefactor in Eq. (^) 
while the short dashed line shows the ln(274)-contribution to 
A. The solid line in (b) represents the analytical result ac- 
cording to Eqs. (^)-(^), while the long dashed line is obtained 
by exact numerical evaluation of Eq. (0) or Eq. (^. The 
short dashed line in (b) shows the ln(2A)^contribution, while 
the dot-dashed line gives the 1/yl '^-contribution. The DM in 
(b) has been calculated on the particular representative node 
characterized by v^^ = Vy = = 9. 



to be estimated. It is clearly seen in Figs. |](c) and §(d) 
that the major A-dependent contribution comes from the 
ln(2^)-term in Eq. (||), while the role of the second term, 
'-^ \/A^, in Eq. (||) is not important at all, especially in 
comparison with the significant A-independent contribu- 
tion from the last two terms in Eq. (||) (see Fig. ||(b)). 

In the case when one of the indices, e.g. ^ 1, but 
the other ones are small and finite (-D > 2), 



Mil 



exp[-i/,ln(2.4)] n^-. (^) , (6) 



where J stands for the Bessel function, the asymptotic 
form for which, with large argument | p^ , must be used in 
Eq. (|6|). Bearing in mind that J^„(i/xM) « {v^/Ay'^/'^ 
for S> 1, it is easy to see that the power-law prefactor, 
(X v~^/'^^ is restored in Eq. (^). The results for this 
particular case are presented in Figs. |l|(b,d) and H(c,f). 
Again, the contribution of the term ln(2v4) in the inverse 
decay length is dominant, and the inverse decay length 
is now larger than for decay along the main diagonal (cf. 
e.g. Figs, ga) and§(c)). 

Finally, if t'x, J^y 2> 1 and is finite (in 3-D), 



2'KV^Vy 



exp 



4 



A 



(7) 



where v± = Vy^:^ Vy and, for both Bessel functions, 
known asymptotic expressions (for the argument in 3i^^ 
and for the argument and index in J^^ ) must be used 
[ p2[ . These asymptotics for both Bessel functions restore 
the v~^/'^ prefactor for the DM in the same manner as for 
Eq. (6). The energy (A)-dependent contribution to A (see 
Fig. 2(e)), in this case, comes mainly from the asymp- 
totic expression for J^^ which is of exponential type if 
v+/{A,Jir^) < 1 (and of cos- type otherwise) p2| . 

Another important point to discuss is the dependence 
of the typical inverse decay length. A, on the energy pa- 
rameters of the model. The only energy-dependent pa- 
rameter in this model is the parameter A which is sim- 
ply related to the spectral characteristics of the model, 
A = A£/4ii2 = DAEopt/iP^E^), where ASb = ADt 
is the band width, AE'opt = Ae is the optical (direct) 
band gap and the ratio of hopping integrals, (3 = ti2/t, 
which can be connected to the thermal gap, AE'th = 
[A£;2p^ + {l3AEbfY/^ - AEb. The spectral parameters 
(Ai?b, Ai?opt and Ai^th) are not independent values be- 
cause the hopping integrals t^^' in the Hamiltonian (|l|) 
can depend on the orbital energies and obtaining such 
a dependence, e.g. ii2(Ae), requires more accurate anal- 
ysis which is not necessary for our approach (that is why 
we plot Xvs A but not the direct gap width Ae in Fig. ||). 
However, the limiting case, Ae — > 0, can be investigated, 
at least numerically (analytical results are valid only for 
A > D) under the reasonable assumption that ti2 ^ 
as Ae ^ 0, and the results are presented in Fig. ||(a), 
confirming the linear scaling of the inverse decay length, 
Aoc^cx A£for^->0 |]. 

This model can be analysed analytically only for a 
s.c. lattice and thus, strictly speaking, we stick to the 
atomic orbital representation. However, numerics show 
that the results (exponential decay of the density ma- 
trix and variation of the decay length with parameters 
of the Hamiltonian) are qualitatively the same for differ- 
ent underlying lattices. This is not surprising because 
the asymptotic behaviour at large distances hardly de- 
pends on the local structural details (the local features 
of the lattice can influence the value of A by a factor of 
~ 2; see Figs. Thus, we believe that our results, at 

least qualitatively, are general. Indeed, the inverse de- 
cay length decreases with decreasing A (see Fig. ||), thus 
reflecting the known delocalization tendency of the DM 
with increasing metallicity (in the bonding/antibonding 
representation A ~ ^^2/^1, being the ratio of the 
bonding-antibonding splitting energy and the metallic 
bandwidth energy). Moreover, the analytical expressions 
for the decay length give correct order-of-magnitude es- 



pressions for the spatial decay of the density matrix de- 
cay in insulators. The results have been obtained for 
a tight-binding Hamiltonian defined on the s.c. lattice 
having two orthogonal bare states on each node. Two 
additional assumptions, namely symmetric bands (equal 
in-band transfer integrals) and relatively weak intraband 
hybridization, allow us to derive exact asymptotic results 
for the density matrix in all three dimensions. The main 
features of the analytic solution are: exponential spa- 
tial decay of the density matrix, anisotropy of the decay 
length, the dependence of the decay length on the energy 
parameters of the model, and the existence of a power- 
law prefactor (cx r~-°/^). All the analytic asymptotic re- 
sults have been supported by precise numerical solutions 
of the problem. 
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